//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Implementacao das funcoes da biblioteca c++ para programa de modelagem de
// curva de luz. Versao C++ da biblioteca roche.h para modelar curvas de luz
// da estrela secundaria em Variaveis Cataclismicas. Esta verso deve ser mais
// generica e permitir outras geometrias de forma mais simples que a versao
// C. Por exeplo, calcular curva de luz para diferentes geometrias, disco,
// estrela esferica, estrela preenchendo lobo de roche e estrela nao 
// preenchendo o lobo de roche. 
//
// Os cabecalhos das funcoes se encontram no arquivo 'light_curve.h'. Com
// ele deve ser possivel de entender o funcionamento da biblioteca sem a 
// nescessidade de ler os codigos deste arquivo.
//
// Definicao da classe Roche para curvas de luz de estrelas preenchendo o 
// Lobo de Roche.
//
// Tiago Ribeiro -- LinuxC++ -- UFSC -- 2006
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+


#include "light_curve.h"
#include "cxx/generic_gcc/include/healpix_base.cc"
#include "cxx/generic_gcc/include/healpix_map.h"
#include "cxx/generic_gcc/include/healpix_map_fitsio.h"
#include "dataFileReader_v0.1.1/dataFileReader.h"
//#include <omp.h>
//#include "lib_zabot/aborta.h"
//#include "erro_handler.h"
  
using namespace std;

// lembrar de incluir light_curve.cpp

Roche::Roche()
{
  
}

Roche::Roche(char *pixfile, const vector<vector<double> > &data_par):Light_Curve(data_par)
{
  roche_par.i = data_par[0][0];
  roche_par.q = data_par[0][2]/data_par[0][1];
  roche_par.f_fac = 1; //para os programas que nao estao atualizados
  roche_par.obs_grav = data_par[1][1];

  //monta  superficie da estrela
  lagrangeana(); // acha rl1, nao trocar ordem!
  roche_par.g_pole = gpole();

  load_pix_file(pixfile);
  faz_super();

	cout << "flux(0) = " << flux(0.0) << endl;

  if(par[1][0] > 0) par[1][0] = 1./flux(0);

  cout << "flux(0) = " << flux(0) << endl;
}

Roche::Roche(const int dec, const vector<vector<double> > &data_par) : Light_Curve(data_par)
{
  // testar se o numero de parametros esta certo!!!

  //  if((data_par[0].size()+data_par[1].size()) != N_PAR_ORBITAL+N_PAR_ATMOSF+int(data_par[5]+data_par[6])-2){
    //    throw Par_Number(data_par.size(), N_PAR_ORBITAL+N_PAR_ATMOSF+int(data_par[5]+data_par[6])-2);
  // cerr << "ERRO: Numero de parametros!\n"
  //	 << "ERRO: Encerrando programa...\n";
  //  exit(0);
  // }

  //reserva espaco para vetores da superficie

  try{
    //int i = 0;
    star.theta.reserve(dec*dec*dec);
    star.phi.reserve(dec*dec*dec);
    star.px.reserve(dec*dec*dec);
    star.py.reserve(dec*dec*dec);
    star.pz.reserve(dec*dec*dec);
    star.nx.reserve(dec*dec*dec);
    star.ny.reserve(dec*dec*dec);
    star.nz.reserve(dec*dec*dec);
    star.I.reserve(dec*dec*dec);
    star.vx.reserve(dec*dec*dec);
    star.vy.reserve(dec*dec*dec);
    star.vz.reserve(dec*dec*dec);
  }
  catch(bad_alloc){
    cout << "Memoria nao pode ser alocada!\n";
    exit(1);
  }

  roche_par.i = data_par[0][0];
  roche_par.q = data_par[0][2]/data_par[0][1];
  if(data_par[1].size() < 6) {
    roche_par.f_fac = 1.; //para os programas que nao estao atualizados e dispara WARNING...
    cout << "WARNING: Programa nao conta com parametro f_fac.\n"
	 << "\tFavor atualizar!\n\tf_fac = 1.\n";
  }
  else{roche_par.f_fac = data_par[1][5];}
  roche_par.obs_grav = data_par[1][1];

  //monta  superficie da estrela
  lagrangeana(); // acha rl1, nao trocar ordem!
  roche_par.g_pole = gpole();
  
  //v_ang(200);

  healpix_pixelization(dec);
  __Dec__=dec;
  faz_super();
  //cout << "# 0 - flux(0) = " << flux(0) << " | par[1][0] = " << par[1][0] << endl;
  
  F20 = par[1][0];

  if(F20 > 0){
    par[1][0] = F20*F20/(flux(0));
  }

  /*  cout << "# 1 - flux(0) = " << flux(0) << " | par[1][0] = " << par[1][0] << endl;
  vector<double> lixo(3);
  lixo[0] = -500e5;
  lixo[1] = 500e5;
  lixo[2] = 50e5;

  cout << rad_prof(0, lixo);
  exit(0);*/

}

Roche::Roche(const surface &StarSurf, const vector<vector<double> > &data_par):Light_Curve(data_par)
{

  /*
  try{
    star.theta (StarSurf.theta);
    star.phi.push_back(StarSurf.phi);
    star.px.push_back(StarSurf.px);
    star.py.push_back(StarSurf.py);
    star.pz.push_back(StarSurf.pz);
    star.nx.push_back(StarSurf.nx);
    star.ny.push_back(StarSurf.ny);
    star.nz.push_back(StarSurf.nz);
    star.I.push_back(StarSurf.I);
    star.vx.push_back(StarSurf.vx);
    star.vy.push_back(StarSurf.vy);
    star.vz.push_back(StarSurf.vz);
  }
  catch(bad_alloc){
    cout << "Memoria nao pode ser alocada!\n";
    exit(1);
  }
  */

  star = StarSurf;

  roche_par.i = data_par[0][0];
  roche_par.q = data_par[0][2]/data_par[0][1];
  if(data_par[1].size() < 6) {
    roche_par.f_fac = 1.; //para os programas que nao estao atualizados e dispara WARNING...
    cout << "WARNING: Programa nao conta com parametro f_fac.\n"
	 << "\tFavor atualizar!\n\tf_fac = 1.\n";
  }
  else{roche_par.f_fac = data_par[1][5];}
  roche_par.obs_grav = data_par[1][1];

  lagrangeana(); // acha rl1, nao trocar ordem!
  roche_par.g_pole = gpole();

  F20 = par[1][0];

  for(int i = 0; i < star.I.size(); i++)
    star.I[i] = 1.0;

  if(F20 > 0){
    par[1][0] = F20*F20/(flux(0));
  }

  for(int i = 0; i < star.I.size(); i++)
    star.I[i] = StarSurf.I[i];


  //if(par[1][0] > 0) par[1][0] = 1./flux(0);

}

double Roche::gpole(void)
{
  double qq, rlob, x, nr, n_r, nx, np;
  
  qq = 1./roche_par.q;
  rlob = ache_r(0, M_PI/2);
  x = rlob;
  n_r = abs(1-x);
  nr = x*x*x;
  n_r = n_r*n_r*n_r;
  nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
  np = abs(nx);
  return np;
}


void Roche::lagrangeana(void)
{ 
  double a=0, b=0, c=0, tfi=0, x=0, x1=0, x2=1, dfi=0,q;
  int i=0, j=0, k=0;

  x1 = 0;
  x2 = 1;
  q = 1./roche_par.q; //para lobo da secundaria

  for(x=0; i<=20 ; i++){
    
      x = ((x2-x1)/2)+x1;
      a = -2*x/((1+q)*pow((x*x),3/2));
      b = -2*q*(x-1)/((1+q)*pow(pow(x-1,2),3/2));
      c = 2*(x-(q/(1+q)));
      dfi = a+b+c;
 
      if(dfi>0){ 
	x2 = x;
      }

      if(dfi<0){
	x1 = x;
      }

      if(dfi>=-1e-5 && dfi<=1e-5){
	rl1 = x;
	break;
      }

    }

}

void Roche::load_pix_file(const char *pixfile)
{
  DataFileReader dfr_pix(pixfile);

  cout << "# Lendo " << pixfile << endl;

  while(!++dfr_pix){
    if(dfr_pix.numberOfColumns() != 2){
      ostringstream erro;
      erro << "(Linha: "<< dfr_pix.lineNumber() << ") - Numero de colunas errado em " 
	   << pixfile << " ...\n"
	   << "Esperando 'theta phi' (2 culunas) - lido " << dfr_pix.numberOfColumns() << endl;
      aborta(erro.str());
    }

	  double theta = dfr_pix.readColumn<double>(1);
	  double phi = dfr_pix.readColumn<double>(2);

    //dfr_pix >> theta >> phi;
    ++dfr_pix;
    theta -= M_PI/2.;

    star.theta.push_back(theta); //theta
    star.phi.push_back(phi); //phi
    star.px.push_back(sin(theta));
    star.py.push_back(cos(theta)*cos(phi));
    star.pz.push_back(cos(theta)*sin(phi));

    cout << sin(theta) << " " 
	 << cos(theta)*cos(phi) << " "
	 << cos(theta)*sin(phi) << endl;

  }
  dfr_pix.close();
}

void Roche::healpix_pixelization(const int dec)
{
  Healpix_Map<double> pixmap(dec, RING);

  pixmap.fill(1);
  
  //cout << "# PIXELIZATION LOG: Npix = " << pixmap.Npix() << endl;

  for(int i = 0; i < pixmap.Npix(); i++){
    pointing vec = pixmap.pix2ang(i);

    //Meu theta eh o phi de Healpix_map e vice versa...
    double theta = vec.phi; //-M_PI/2.0;
    double phi = vec.theta-M_PI/2.0;
    if(star.theta.size() < pixmap.Npix()){

#ifdef DEBUG
      if(i == 0) cout << "# Incrementing star.theta structure\n";
#endif
      star.theta.push_back(theta); //theta
      star.phi.push_back(phi); //phi
      star.px.push_back(sin(phi));
      star.py.push_back(cos(phi)*cos(theta));
      star.pz.push_back(cos(phi)*sin(theta));
      star.nx.push_back(0.0);
      star.ny.push_back(0.0);
      star.nz.push_back(0.0);
      star.vx.push_back(0.0);
      star.vy.push_back(0.0);
      star.vz.push_back(0.0);
    } 
    else{
      star.theta[i]=(theta); //theta
      star.phi[i]=(phi); //phi
      star.px[i]=(sin(phi));
      star.py[i]=(cos(phi)*cos(theta));
      star.pz[i]=(cos(phi)*sin(theta));
    }
  }
  star.I.resize(star.theta.size());
_size_=star.theta.size();
}

void Roche::healpixmap_2_vtk(int dec, ostringstream &vtk_out)
{
	cout << "# UNDEFINED\n";
/*
	//colocando cabecalhos...
  vtk_out << "# vtk DataFile Version 2.0\n"
	  << "Data test\n"
	  << "ASCII\n"
	  << "DATASET UNSTRUCTURED_GRID\n"
	  << "POINTS ";

	//ealpix_Map<double> pixmap(__Dec__,RING);
	//pixmap.fill(1);
	//if(pixmap.Npix != star.theta.size())
	//	aborta("HEALPIXMAP_2_VTK(): Superficie definida de maneira diferente do esperado");

	//Numero de pixel 3 vezes maior q o numero de pixels de verdade... 
	//Pegando as bordas de cada pixel
  
  surface pixmap;
  openGL_image(dec,pixmap);
  
  vtk_out << pixmap.theta.size()*3 << " float\n";
	//Escrevo a posicao de cada borda do pixel em ordem

	//ESCREVER CODIGO!!!

	//Definicao das celulas... 1 celula para cada pixel!
	// CELLS 'numero de pixels' 'numero de parametros para ler (4 por pixel)'
  vtk_out << "\nCELLS " << pixmap.NPIX() << " " << pixmap.NPIX()*4 << endl;
	
  int cont0=0;
  for(int i = 0; i < pixmap.NPIX()*3; i++){
		//vtk_out << 
		cont0++;
		if(cont0==3) {
			cont0=0;
			vtk_out<<endl;
		}
	}

	vtk_out << endl << "CELLS " << pixmap.Npix() << " " pixmap.Npix()*4 << endl;
	
	for(int i = 0; i < pixmap.Npix(); i+=3)
		vtk_out << "3 " << i+0 << " " << i+1 << " " << i+2 << endl;
	
	vtk_out << "CELL_TYPES " << pixmap.Npix() << endl;

	for(int i = 0; i < pixmap.Npix(); i++) vtk_out << "5" << endl;
	vtk_out<<endl<<endl;

	vtk_out << "POINT_DATA " << pixmap.Npix*3 << endl
		<< "SCALARS intensity float 1\n"
		<< "LOOKUP_TABLE default\n";

	for(int i = 0; i < pixmap.Npix(); i++)
		vtk_out << star.I[i] << " "  << star.I[i] << " " << star.I[i] << endl;
	
*/
}

void Roche::v_ang(const int dec)
{
  double v, u, a1, dvv, AZIM, vvmin, vvmax, rlk;
  const double vmax = M_PI/2.;
  const double vmin = -M_PI/2.;
  const double dv = M_PI/double(dec);
  const double vmax2 = (M_PI-dv)/2.;
  const double vmin2 = -vmax2;
  const double rl0 = ache_r(0,0);
  int i=0, j=0, k=0,k2;

  // mudanca de definicao para argumento da funcao! LEMBRE-SE
  AZIM = double(dec); 
  v = vmax;

  //coloca os dois polos
  v = 2.*M_PI;
  u = vmax;
 
  star.theta.push_back(v); //theta
  star.phi.push_back(u); //phi
  a1 = cos(u);
  star.px.push_back(sin(u));
  star.py.push_back(a1*cos(v));
  star.pz.push_back(a1*sin(v));

  u = vmin;
  star.theta.push_back(v); //theta
  star.phi.push_back(u); //phi
  a1 = cos(u);
  star.px.push_back(sin(u));
  star.py.push_back(a1*cos(v));
  star.pz.push_back(a1*sin(v));
  

  //nao repetir os polos 
  //vmin+=(dv/2.);
  //vmax-=(dv/2.);

  try{
    for(u = vmax2; u >= vmin2; u -= dv ){
      
      if (u == vmax || u == vmin) i = 1;
      else {
	dvv = ceil(AZIM*cos(u)*rl0/ache_r(0,u));
	dvv = M_PI/dvv;
      }
      
      a1 = cos(u);
      
      vvmin = -M_PI/2.+dvv/2.;
      vvmax = M_PI/2.-dvv/2.; //-dvv/2.;
      
      for(v = vvmin; v < vvmax; v+=dvv){
	
	star.theta.push_back(v); //theta
	star.phi.push_back(u); //phi
	star.px.push_back(sin(u));
	star.py.push_back(a1*cos(v));
	star.pz.push_back(a1*sin(v));
	
	star.theta.push_back(M_PI+v); //theta
	star.phi.push_back(u); //phi
	star.px.push_back(sin(u));
	star.py.push_back(a1*cos(M_PI+v));
	star.pz.push_back(a1*sin(M_PI+v));

      }
      // coloco v = vvmax = M_PI/2.-dvv/2. a mao...


    }
  }
  catch (bad_alloc){ 
    cerr << "Mem nao pode ser alocada!\n";
    exit(0);
  } 

  cout << "# PIXELIZATION LOG: Npix = " << star.theta.size() << endl;
}

double Roche::ache_r(double v0, double u0)
{
  double q, fi, fimax, R=0, R1, R2, w, a, b, c, d, r_2, x, y;
  int i=0;

  q = 1./roche_par.q; // lobo da secundaria

  R1 = 0.;
  R2 = rl1;	
  a = 2./((1+q)*R2);
  b = 2.*q/((1+q)*sqrt(pow(R2,2.)-(2.*R2)+1));
  c = pow((R2)-(q/(1+q)),2);
  /*d = pow(R2,2);*/ 
  fimax = (a+b+c)*roche_par.f_fac;

  for(R=0 ; R1 < R2 ;){
    R = ((R2-R1)/2)+R1; 
    r_2 = sqrt(pow(R,2.)-2.*R*sin(u0)+1);
    x = R*sin(u0);
    y = R*cos(u0)*cos(v0);	
    a = 2./((1+q)*R);
    b = 2.*q/((1+q)*r_2);
    c = pow(x-(q/(1+q)),2.);
    d = pow(y,2.);
    fi = a+b+c+d;
   
    
    if(fi < fimax){ 
      R2 = R;
    }
    if(fi > fimax){
      R1 = R;
    }
    if(fabs(fi-fimax) <= 1e-6){
      /*fprintf(vp," %e %e %e %e %e\n", fimax, fi, R, R1, R2);*/ 
      return R; 
    }
    
  }

  cout << "O ponto " << fi << ' ' << R << ' ' << u0 << ' ' << v0 
       << " nao foi encontrado!\nO programa sera terminado!\n";
  exit(1);

}

double Roche::ache_r2(double v0, double u0)
{
  double q, fi, fimax, R=0, R1, R2, w, a, b, c, d, r_2, x, y;
  int i=0;

  q = 1./roche_par.q; // lobo da secundaria

  R1 = 0.;
  R2 = rl1;	
  a = 2./((1+q)*R2);
  b = 2.*q/((1+q)*sqrt(pow(R2,2.)-(2.*R2)+1));
  c = pow((R2)-(q/(1+q)),2);
  /*d = pow(R2,2);*/ 
  fimax = (a+b+c)*roche_par.f_fac;

  for(R=0 ; R1 < R2 ;){
    R = ((R2-R1)/2)+R1; 
    x = R*cos(u0)*sin(v0); //sin(u0);
    y = R*cos(u0)*cos(v0);
    r_2 = sqrt(pow(R,2.)-2.*x+1);	
    a = 2./((1+q)*R);
    b = 2.*q/((1+q)*r_2);
    c = pow(x-(q/(1+q)),2.);
    d = pow(y,2.);
    fi = a+b+c+d;
   
    
    if(fi < fimax){ 
      R2 = R;
    }
    if(fi > fimax){
      R1 = R;
    }
    if(fabs(fi-fimax) <= 1e-6){
      /*fprintf(vp," %e %e %e %e %e\n", fimax, fi, R, R1, R2);*/ 
      return R; 
    }
    
  }

  cout << "O ponto " << fi << ' ' << R << ' ' << u0 << ' ' << v0 
       << " nao foi encontrado!\nO programa sera terminado!\n";
  exit(1);

}

void Roche::CalcPixSize(void){
	// Calculando o raio da estrela 
	  using namespace Fis_const;
	double q = roche_par.q;
	double ww = 1./(par[0][3]*3600.); //frequencia angular
	//sep. orbital
	double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4.0*M_PI*M_PI*ww*ww), 1./3.);
	
	// Raio segundo formula de Eggleton (1983)
	double rl2 = 0.49*pow(q,2./3.)/(0.69*pow(q,2./3.)+log(1+pow(q, 1./3.))) * a * 2.*M_PI*ww; //6.955e10;
	// Area total da secundaria
	double At2 = 2.*M_PI*rl2*rl2/3.;
	// Area do pixel = At2 / npix
	pix_size(At2/npix()); //sqrt()*2.*M_PI*ww*a*rl2/1e5
}


void Roche::faz_super(void)
{
  int k = 0;
  double u, u1, v, rlob,pu, pb, rm, qq, If, gp;
  double x, y, z, nr, n_r, nx, ny, nz, np, n_s, lt_o=0, betha;
  double *fs;

  using namespace Fis_const;
  //testando superficie da estrela
  if((star.nx.size() != star.ny.size()) || (star.nz.size() != star.ny.size()) || (star.vx.size() != star.ny.size()) || (star.vy.size() != star.vz.size()) || (star.ny.size() != star.vz.size()))
    aborta("-> Problema com a definicao da superficie da estrela");
  //funcao que cria a superficie da estrela. Ela guarda raio, versores e
  //brilho de TODOS ponto da superficie.

  lt_o = 0;
  qq = 1./roche_par.q; // lobo da secundaria
  k = 0;

  double ww = 1./(par[0][3]*3600.); //frequencia angular
  //sep. orbital
  double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4.0*M_PI*M_PI*ww*ww), 1./3.);
  double r2 = 1./(par[0][2]/par[0][1]+1.);//raio orbital da sec. em unidade de a
// cout << "# Velocidade Orbital da secundaria= " << 2.*M_PI*ww*a*r2/1e5*sin(par[0][0]) <<endl;
 // cout << "# Raio orbital da secundaria (em a) = " << r2 << endl;
 // cout << "# Raio da estrela (em Rsun) = " << ache_r(0., 0.)*a/6.955e10 << endl;
  for(int w = 0; w < star.theta.size(); w++) {
    rlob = ache_r(star.theta[w], star.phi[w]);
    x = star.px[w]*rlob;
    y = star.py[w]*rlob;
    z = star.pz[w]*rlob;
    nr = sqrt(x*x+y*y+z*z);
    n_r = sqrt((1-x)*(1-x)+(y*y)+(z*z));
    nr = nr*nr*nr;
    n_r = n_r*n_r*n_r;
    nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
    ny = -y*(1./nr+qq/n_r-(qq+1));
    nz = -z*(1./nr+qq/n_r);
    np = sqrt(nx*nx+ny*ny+nz*nz);

    //    if( fabs(x) == 1 ) betha = x > 0 ? 0. : M_PI;
    //else betha = atan(y/x); 
    //pu = (*(p[1]+w)-M_PI/2.)/(M_PI/6.);
    //pb = (betha-d[3])/(M_PI);
    
    /* 
       !SPOT MAKER! -> colocar numa funcao separada!
    pu = (M_PI/2-(star.phi[w]))-d[3])/(M_PI/6.);
    if(star.theta[w] <= M_PI) pb = star.theta[w]*d[3]/(M_PI/6.);
    else pb = (2*M_PI- star.theta[w])*d[3]/(M_PI/6.);
    Is = pow(np/o[4], o[2])*(1+(d[2]*exp(-(pu*pu+pb*pb))));
       !SPOT MAKER! 
    */

    //Intesidade do pixel atual. Com obscurecimento gravitacional
    // o[2] -> coeficiente de obscurecimento gravitacional.

    star.I[w] = pow(np/roche_par.g_pole, roche_par.obs_grav);

    star.px[w]=x; //rlob;
    star.py[w]=y; //rlob;
    star.pz[w]=z; //rlob;

    //    cout << star.theta[w] << " "
    //	 << star.phi[w] << " "
    //	 << star.px[w] << " "
    //   << star.py[w] << " "
    //	 << star.pz[w] << endl;
 
    //Armazenando tudo!!
    //    printf("%lf %lf\n",np,o[4]);
    //if(w >= star.nx.size()){
    star.nx[w]=(nx);
    star.ny[w]=(ny);
    star.nz[w]=(nz);
    star.vx[w]=(-2.*ww*y*M_PI*a); //em km/s
    star.vy[w]=(2.*ww*x*M_PI*a); //em km/s
    if(vorb()) star.vy[w]+=-2.*M_PI*ww*a*r2; //adiciona velocidade orbital
    star.vz[w]=(0.);
    //}

    //printf("-->%f %f %f %f %f\n",*(fs-5), *(fs-4),*(fs-3),*(fs-2),*(fs-1));
    //printf("<--%f %f %f %f %f\n",rlob, nx, ny, nz, If);
    //Obscurecimento de borda. Calculado na hora!
    //o[0] -> coeficiente de obscurecimento de borda.
    //o[1] -> coeficiente nao linear de obscurecimento de borda.
    //If = If*(1-o[0]*(1-n_s)-o[1]*(1-sqrt(n_s)));
  }
	
  
}

void Roche::refaz_super(void)
{
  int k = 0;
  double u, u1, v, rlob,pu, pb, rm, qq, If, gp;
  double x, y, z, nr, n_r, nx, ny, nz, np, n_s, lt_o=0, betha;
  double *fs;

  using namespace Fis_const;
  //testando superficie da estrela
  if((star.nx.size() != star.ny.size()) || (star.nz.size() != star.ny.size()) || (star.vx.size() != star.ny.size()) || (star.vy.size() != star.vz.size()) || (star.ny.size() != star.vz.size()))
    aborta("-> Problema com a definicao da superficie da estrela");
  //funcao que cria a superficie da estrela. Ela guarda raio, versores e
  //brilho de TODOS ponto da superficie.

  lt_o = 0;
  qq = 1./roche_par.q; // lobo da secundaria
  k = 0;

  double ww = 1./(par[0][3]*3600.); //frequencia angular
  //sep. orbital
  double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4*M_PI*M_PI*ww*ww), 1./3.);
  double r2 = 1./(par[0][2]/par[0][1]+1.);//raio orbital da sec. em unidade de a

  for(int w = 0; w < star.theta.size(); w++) {
    rlob = ache_r(star.theta[w], star.phi[w]);
    x = star.px[w];// *rlob;
    y = star.py[w];// *rlob;
    z = star.pz[w];// *rlob;
    nr = sqrt(x*x+y*y+z*z);
    n_r = sqrt((1-x)*(1-x)+(y*y)+(z*z));
    nr = nr*nr*nr;
    n_r = n_r*n_r*n_r;
    nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
    ny = -y*(1./nr+qq/n_r-(qq+1));
    nz = -z*(1./nr+qq/n_r);
    np = sqrt(nx*nx+ny*ny+nz*nz);

    //    if( fabs(x) == 1 ) betha = x > 0 ? 0. : M_PI;
    //else betha = atan(y/x); 
    //pu = (*(p[1]+w)-M_PI/2.)/(M_PI/6.);
    //pb = (betha-d[3])/(M_PI);
    
    /* 
       !SPOT MAKER! -> colocar numa funcao separada!
    pu = (M_PI/2-(star.phi[w]))-d[3])/(M_PI/6.);
    if(star.theta[w] <= M_PI) pb = star.theta[w]*d[3]/(M_PI/6.);
    else pb = (2*M_PI- star.theta[w])*d[3]/(M_PI/6.);
    Is = pow(np/o[4], o[2])*(1+(d[2]*exp(-(pu*pu+pb*pb))));
       !SPOT MAKER! 
    */

    //Intesidade do pixel atual. Com obscurecimento gravitacional
    // o[2] -> coeficiente de obscurecimento gravitacional.

    star.I[w] = pow(np/roche_par.g_pole, roche_par.obs_grav);

    star.px[w]=x; //rlob;
    star.py[w]=y; //rlob;
    star.pz[w]=z; //rlob;

    //    cout << star.theta[w] << " "
    //	 << star.phi[w] << " "
    //	 << star.px[w] << " "
    //   << star.py[w] << " "
    //	 << star.pz[w] << endl;
 
    //Armazenando tudo!!
    //    printf("%lf %lf\n",np,o[4]);
    //if(w >= star.nx.size()){
    star.nx[w]=(nx);
    star.ny[w]=(ny);
    star.nz[w]=(nz);
    star.vx[w]=(2.*ww*y*M_PI*a); //em km/s
    star.vy[w]=(-2.*ww*x*M_PI*a); //em km/s
    if(vorb()) star.vy[w]+=2.*M_PI*ww*a*r2; //adiciona velocidade orbital
    star.vz[w]=(0.);
    //}

    //printf("-->%f %f %f %f %f\n",*(fs-5), *(fs-4),*(fs-3),*(fs-2),*(fs-1));
    //printf("<--%f %f %f %f %f\n",rlob, nx, ny, nz, If);
    //Obscurecimento de borda. Calculado na hora!
    //o[0] -> coeficiente de obscurecimento de borda.
    //o[1] -> coeficiente nao linear de obscurecimento de borda.
    //If = If*(1-o[0]*(1-n_s)-o[1]*(1-sqrt(n_s)));
  }
  
}


void Roche::openGL_Image(const int dec, surface &Image)
{
  //pixerizao para OpenGl - npix constante para todas as dec
  GL_pix(dec, Image);

  //determinando superficie de roche...

  faz_super(Image);

/*  double fase = 0.25;
  double u1 = -fase*2.*M_PI;
  
  double z1 = cos(par[0][0]);
  double z2 = sin(par[0][0]);
  double f1 = sin(u1);
  double f2 = cos(u1);
  
  double x1 = z2*f2; //sin(par[0][0])*cos(u1);
  double y1 = f1*z2; //sin(u1)*sin(par[0][0]);
  int nobs = int(par[1][2]);
*/
  //vector<int> Itrac(Image.theta.size(), 1);

  //Convertendo imagem real para roche...
  //cout << "aqui 0\n";
  Image.I.resize(Image.px.size());
//#pragma omp parallel for
  for(int i = 0; i < Image.px.size(); i++){
    int index = proc_pix(Image.px[i], Image.py[i], Image.pz[i]);
/*    double n_s = star.nx[index]*x1+star.ny[index]*y1+star.nz[index]*z1;
    n_s = n_s/sqrt(star.nx[index]*star.nx[index]+star.ny[index]*star.ny[index]+star.nz[index]*star.nz[index]);
*/
    Image.I[i] = star.I[index]; //*obscurecimento_de_borda(fabs(n_s),nobs);;
    //cout << Image.px[i] << " " << Image.py[i] << " " << Image.pz[i] << endl;
    //Itrac[i]++;
  }
  //cout << "aqui 1\n";
  //cout << endl;
  //faz_super2(Image);
  //for(int i = 0; i < star.theta.size(); i++)    
  //  cout << star.px[i] << " " << star.py[i] << " " << star.pz[i] << endl;
  //for(int i = 0; i < Image.theta.size(); i++)
  //  Image.I[i]/=double(Itrac[i]);
}

void Roche::openGL_Image_polar(const int dec, surface &Image)
{
  //pixerizao para OpenGl - npix constante para todas as dec
  GL_pix_polar(dec, Image);

  //determinando superficie de roche...

  //faz_super2(Image);

  //vector<int> Itrac(Image.theta.size(), 1);

  //Convertendo imagem real para roche...
  //cout << "aqui 0\n";
  Image.I.resize(Image.px.size());
  for(int i = 0; i < Image.px.size(); i++){
    int index = proc_pix(Image.px[i], Image.py[i], Image.pz[i]);
    Image.I[i] = star.I[index];
    //cout << Image.px[i] << " " << Image.py[i] << " " << Image.pz[i] << endl;
    //Itrac[i]++;
  }
  //cout << "aqui 1\n";
  //cout << endl;
  //faz_super2(Image);
  //for(int i = 0; i < star.theta.size(); i++)    
  //  cout << star.px[i] << " " << star.py[i] << " " << star.pz[i] << endl;
  //for(int i = 0; i < Image.theta.size(); i++)
  //  Image.I[i]/=double(Itrac[i]);
}

void Roche::GL_pix(const int dec, surface &Image)
{
  const double dphi = M_PI/double(dec);
  const double dtheta = 2*M_PI/double(dec);

  double phi1;
  double phi2;

  for(int j = 0; j < dec; j++){
    phi1 = -M_PI/2+double(j)*dphi;
    phi2 = -M_PI/2+double(j+1)*dphi;
    //phi1 = double(j)*dphi;
    //phi2 = double(j+1)*dphi;
    double theta;
    for(int i = 0; i <= dec; i++){

      theta = double(i)*dtheta;
 
      Image.theta.push_back(theta);
      Image.phi.push_back(phi1);
      Image.px.push_back(sin(phi1));
      Image.py.push_back(cos(phi1)*cos(theta));
      Image.pz.push_back(cos(phi1)*sin(theta));
     
      Image.theta.push_back(theta);
      Image.phi.push_back(phi2);
      Image.px.push_back(sin(phi2));
      Image.py.push_back(cos(phi2)*cos(theta));
      Image.pz.push_back(cos(phi2)*sin(theta));
      //cout << sin(phi2) << " " 
      //     << cos(phi2)*cos(theta) << " "
      //     << cos(phi2)*sin(theta) << endl;
    }
  }


}

void Roche::GL_pix_polar(const int dec, surface &Image)
{
  const double dphi = M_PI/double(dec);
  const double dtheta = 2*M_PI/double(dec);

  double phi1;
  double phi2;

  for(int j = 0; j < dec; j++){
    phi1 = -M_PI/2+double(j)*dphi;
    phi2 = -M_PI/2+double(j+1)*dphi;
    //phi1 = double(j)*dphi;
    //phi2 = double(j+1)*dphi;
    double theta;
    for(int i = 0; i <= dec; i++){

      theta = double(i)*dtheta;
 
      Image.theta.push_back(theta);
      Image.phi.push_back(phi1);
      Image.pz.push_back(sin(phi1));
      Image.py.push_back(cos(phi1)*cos(theta));
      Image.px.push_back(cos(phi1)*sin(theta));
     
      Image.theta.push_back(theta);
      Image.phi.push_back(phi2);
      Image.pz.push_back(sin(phi2));
      Image.py.push_back(cos(phi2)*cos(theta));
      Image.px.push_back(cos(phi2)*sin(theta));
      //cout << sin(phi2) << " " 
      //     << cos(phi2)*cos(theta) << " "
      //     << cos(phi2)*sin(theta) << endl;
    }
  }


}


void Roche::faz_super(surface &Image)
{
  int k = 0;
  double u, u1, v, rlob,pu, pb, rm, qq, If, gp;
  double x, y, z, nr, n_r, nx, ny, nz, np, n_s, lt_o=0, betha;
  double *fs;

  using namespace Fis_const;

  //funcao que cria a superficie da estrela. Ela guarda raio, versores e
  //brilho de TODOS ponto da superficie.

  lt_o = 0;
  qq = 1./roche_par.q; // lobo da secundaria
  k = 0;

  double ww = 1./(par[0][3]*3600.); //frequencia angular
  //sep. orbital
  double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4*M_PI*M_PI*ww*ww), 1./3.);
  double r2 = 1./(par[0][2]/par[0][1]+1.);//raio orbital da sec. em unidade de a

  for(int w = 0; w < Image.theta.size(); w++) {
    rlob = ache_r(Image.theta[w], Image.phi[w]);
    x = Image.px[w]*rlob;
    y = Image.py[w]*rlob;
    z = Image.pz[w]*rlob;
    if(w == 0) {
      x = rlob;
      Image.px[w] = rlob;
    }
    nr = sqrt(x*x+y*y+z*z);
    n_r = sqrt((1-x)*(1-x)+(y*y)+(z*z));
    nr = nr*nr*nr;
    n_r = n_r*n_r*n_r;
    nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
    ny = -y*(1./nr+qq/n_r-(qq+1));
    nz = -z*(1./nr+qq/n_r);
    np = sqrt(nx*nx+ny*ny+nz*nz);

    //    if( fabs(x) == 1 ) betha = x > 0 ? 0. : M_PI;
    //else betha = atan(y/x); 
    //pu = (*(p[1]+w)-M_PI/2.)/(M_PI/6.);
    //pb = (betha-d[3])/(M_PI);
    
    /* 
       !SPOT MAKER! . colocar numa funcao separada!
    pu = (M_PI/2-(Image.phi[w]))-d[3])/(M_PI/6.);
    if(Image.theta[w] <= M_PI) pb = Image.theta[w]*d[3]/(M_PI/6.);
    else pb = (2*M_PI- Image.theta[w])*d[3]/(M_PI/6.);
    Is = pow(np/o[4], o[2])*(1+(d[2]*exp(-(pu*pu+pb*pb))));
       !SPOT MAKER! 
    */

    //Intesidade do pixel atual. Com obscurecimento gravitacional
    // o[2] . coeficiente de obscurecimento gravitacional.

    //Image.I.push_back(pow(np/roche_par.g_pole, roche_par.obs_grav));
    Image.I.push_back(1);
 
    Image.px[w]*=rlob;
    Image.py[w]*=rlob;
    Image.pz[w]*=rlob;

    //Armazenando tudo!!
    //    printf("%lf %lf\n",np,o[4]);
    Image.nx.push_back(nx);
    Image.ny.push_back(ny);
    Image.nz.push_back(nz);

    Image.vx.push_back(-2.*ww*y*M_PI*a);
    Image.vy.push_back(2.*ww*x*M_PI*a);

    //printf("-.%f %f %f %f %f\n",*(fs-5), *(fs-4),*(fs-3),*(fs-2),*(fs-1));
    //printf("<--%f %f %f %f %f\n",rlob, nx, ny, nz, If);
    //Obscurecimento de borda. Calculado na hora!
    //o[0] . coeficiente de obscurecimento de borda.
    //o[1] . coeficiente nao linear de obscurecimento de borda.
    //If = If*(1-o[0]*(1-n_s)-o[1]*(1-sqrt(n_s)));
  }
  
}

double Roche::R2_rsun(void){
	using namespace Fis_const;

	double qq = 1./roche_par.q; // lobo da secundaria
	
	double ww = 1./(par[0][3]*3600.); //frequencia angular
	//sep. orbital
	double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4.0*M_PI*M_PI*ww*ww), 1./3.);
	
	return ache_r(0., 0.)*a/6.955e10;
	
}

void Roche::faz_super2(surface &Image)
{
  int k = 0;
  double u, u1, v, rlob,pu, pb, rm, qq, If, gp;
  double x, y, z, nr, n_r, nx, ny, nz, np, n_s, lt_o=0, betha;
  double *fs;

  using namespace Fis_const;

  //funcao que cria a superficie da estrela. Ela guarda raio, versores e
  //brilho de TODOS ponto da superficie.

  lt_o = 0;
  qq = 1./roche_par.q; // lobo da secundaria
  k = 0;

  double ww = 1./(par[0][3]*3600.); //frequencia angular
  //sep. orbital
  double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4*M_PI*M_PI*ww*ww), 1./3.);
  double r2 = 1./(par[0][2]/par[0][1]+1.);//raio orbital da sec. em unidade de a

  for(int w = 0; w < Image.theta.size(); w++) {
    rlob = ache_r2(Image.theta[w], Image.phi[w]);
    x = Image.px[w]*rlob;
    y = Image.py[w]*rlob;
    z = Image.pz[w]*rlob;
    
    nr = sqrt(x*x+y*y+z*z);
    n_r = sqrt((1-x)*(1-x)+(y*y)+(z*z));
    nr = nr*nr*nr;
    n_r = n_r*n_r*n_r;
    nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
    ny = -y*(1./nr+qq/n_r-(qq+1));
    nz = -z*(1./nr+qq/n_r);
    np = sqrt(nx*nx+ny*ny+nz*nz);

    //    if( fabs(x) == 1 ) betha = x > 0 ? 0. : M_PI;
    //else betha = atan(y/x); 
    //pu = (*(p[1]+w)-M_PI/2.)/(M_PI/6.);
    //pb = (betha-d[3])/(M_PI);
    
    /* 
       !SPOT MAKER! . colocar numa funcao separada!
    pu = (M_PI/2-(Image.phi[w]))-d[3])/(M_PI/6.);
    if(Image.theta[w] <= M_PI) pb = Image.theta[w]*d[3]/(M_PI/6.);
    else pb = (2*M_PI- Image.theta[w])*d[3]/(M_PI/6.);
    Is = pow(np/o[4], o[2])*(1+(d[2]*exp(-(pu*pu+pb*pb))));
       !SPOT MAKER! 
    */

    //Intesidade do pixel atual. Com obscurecimento gravitacional
    // o[2] . coeficiente de obscurecimento gravitacional.

    //Image.I.push_back(pow(np/roche_par.g_pole, roche_par.obs_grav));
    Image.I.push_back(1);
 
    Image.px[w]*=rlob;
    Image.py[w]*=rlob;
    Image.pz[w]*=rlob;

    //Armazenando tudo!!
    //    printf("%lf %lf\n",np,o[4]);
    Image.nx.push_back(nx);
    Image.ny.push_back(ny);
    Image.nz.push_back(nz);

    Image.vx.push_back(-2.*ww*y*M_PI*a);
    Image.vy.push_back(2.*ww*x*M_PI*a);
    
    Image.theta[w] = 0;
    Image.phi[w] = 0;
    //printf("-.%f %f %f %f %f\n",*(fs-5), *(fs-4),*(fs-3),*(fs-2),*(fs-1));
    //printf("<--%f %f %f %f %f\n",rlob, nx, ny, nz, If);
    //Obscurecimento de borda. Calculado na hora!
    //o[0] . coeficiente de obscurecimento de borda.
    //o[1] . coeficiente nao linear de obscurecimento de borda.
    //If = If*(1-o[0]*(1-n_s)-o[1]*(1-sqrt(n_s)));
  }
  
}

void Roche::NewI(const gsl_vector &newi)
{
  if(_size_ != newi.size){
    cout << "# - Roche::NewI(const gsl_vector &)\n"
	 << "# -  Operacao nao pode ser realizada, imagens de tamanhos diferentes\n"
	 << "# - star.I[" << star.I.size() << "] | New.I[" << newi.size << "]" << endl;
    aborta("FATAL ERROR");
  }

  for(int i = 0; i < newi.size; i++){
    double &nx = star.nx[i];
    double &ny = star.ny[i];
    double &nz = star.nz[i];
    double np = sqrt(nx*nx+ny*ny+nz*nz);
    star.I[i] = gsl_vector_get(&newi,i)*pow(np/roche_par.g_pole, roche_par.obs_grav);
  }
}

/*
double inline Roche::oculta(vector<double>& p,double fase)
{//Calcula se lobo de roche oculta pixel na posicao x,y,z no vetor pos
  double u1 = -fase*2.*M_PI;
  double z1 = cos(par[0][0]);
  double z2 = sin(par[0][0]);
  double f1 = sin(u1);
  double f2 = cos(u1);
  
  double x1 = z2*f2; //sin(par[0][0])*cos(u1);
  double y1 = z2*f1; //sin(u1)*sin(par[0][0]);

  double x = p[0]*x1;
  double y = y1*p[0]/x1;
  double z = z1*p[0]/x1;
  double r = sqrt(x*x+y*y+z*z);
  //if (r > get_RL1()) return 0.0;
  if(Phi(x,y,z) > Phi0()) {
    return 1.0;
  }
  else {
    return 0.0;
  }

  //return 0;
}
*/

/*
double inline Roche::oculta(vector<double>& p,vector<double> &opos)
{//Calcula se lobo de roche oculta pixel na posicao x,y,z no vetor pos
	
	// rotacionando ponto p em fase e em inclinacao
	double x = p[0]-1.0;
	double y = p[1];
	double z = p[2];
	double yp = y*opos[3]-x*opos[2];
	double xp = x*opos[3]+y*opos[2];
	double zp = z*opos[1]-xp*opos[0];

	if (xp >0 && sqrt(yp*yp + zp*zp) < ache_r(M_PI/2., 0.0)) {
		return 0;
	}
	
	x = p[3];
	y = p[1];
	z = p[2];
	yp = y*opos[3]-x*opos[2];
	xp = x*opos[3]+y*opos[2];
	zp = z*opos[1]-xp*opos[0];
	//cout << yp << " " << zp << " " << xp << endl;	

	
	return 1.0;
}*/
